MONDian three-body predictions for LISA Pathfinder 



O 

(N 

< 
(N 

o 

u 

43 

Or 

i 

o 

CO 

03 



> 
O 

o 

On 
O 



13 



Neil Bevis^'B Joao Magueijo/'Q Christian Trenkel, 2 '[l and Steve Kemble 2 

'Theoretical Physics, Blackett Laboratory, Imperial College, London, SW7 2BZ, United Kingdom 
Astrium Ltd, Gunnels Wood Road, Stevenage SGI 2AS, United Kingdom 
(Dated: August 27, 2010) 

In previous work it was shown that MOND theories predict anomalously strong tidal stresses near 
the saddle points of the Newtonian gravitational potential. An analytical examination of the saddle 
between two bodies revealed a linear and a non-linear solution, valid for the outer and inner regions. 
Here we present a numerical algorithm for solving the MOND equations. We check the code against 
the two-body analytical solutions and explore the region transitioning between them. We then 
develop a a realistic model for the MONDian effects on the saddles of the Sun-Earth-Moon system 
(including further sources is straightforward). For the Sun-Earth saddle we find that the two-body 
results are almost unchanged, with corrections increasing from full to new Moon. In contrast, the 
Moon saddle is an intrinsically three-body problem, but we numerically find a recipe for adapting the 
two-body solution to this case, by means of a suitable re-scaling and axis re-orientation. We explore 
possible experimental scenarios for LISA Pathfinder, and the prospect of a visit to the saddle(s) 
at the end of the mission. Given the chaotic nature of the orbits, awareness of the full range of 
the possibilities is crucial for a realistic prediction. We conclude that even with very conservative 
assumptions on the impact parameter, the accelerometers are abundantly sensitive to vindicate or 
rule out the theory. 



I. INTRODUCTION 

Modified Newtonian dynamics (MOND [l|) is a scheme 
that was first proposed for explaining observed dynami- 
cal properties of galaxies without invoking dark matter. 
The scheme was later incorporated into a proper the- 
ory with a Lagrangian formulation Q, but valid only in 
the non-relativistic regime. Still later a fully covariant 
gravitational theory was found containing MOND phe- 
nomenology as a non-relativistic limit. This theory was 
named TeVeS [|| and alternatives have been proposed 
(e.g. 0-|6j]). The observational features to be studied in 
this paper depend only on their (shared)non-relativistic 
limit. Constraints arising from lensing 0-13' or cosmol- 
ogy [l(| have no bearing here and indicate issues arising 
from the relativistic extension of these theories. 

The opposition between MOND and dark matter 
leaves considerable doubts as to how to interpret new 
astrophysical and cosmological data. A fair comparison 
requires re-evaluating, within each approach, the whole 
set of assumptions underlying the new observations (see 
for example the controversy surrounding the bullet clus- 
ter [llHlJl). For this reason the debate would benefit 
from a direct probe, in the form of a laboratory or So- 
lar System experiment. Such a perspective motivates 
widespread dark matter searches. The analogous "back- 
yard" tests of MOND include searching for anomalies 
in planetary and spacecraft trajectories [L5l - [l7| , stronger 
tidal stresses in the vicinity of saddle points of the New- 
tonian potential (l8j or Solar System manifestations of 



the MOND external field effect [19(. 

In this paper we focus on the MONDian effects on the 
saddle points of the gravitational field, with particular 
emphasis on the region where the Earth and Sun pulls 
cancel (which, we stress, is not at the Lagrange point 
LI). We provide realistic predictions for the future LISA 
Pathfinder spacecraft [2(| , which could plausibly be redi- 
rected to the Earth-Sun saddle, once its primary goals at 
LI have been completed. We also consider the possible 
benefits of re-direction to the saddle point near the Moon 
(there are not two separate Moon-Earth and Moon-Sun 
saddle points) and discuss its related practical issues. 

As already mentioned, the predictions we calculate 
stem from the non-relativistic limit of TeVeS For 
this theory gravity is described by the total potential 
<f> = <f> N + <fi, where $ N is the Newtonian potential and <j> 
is an additional MOND component. The latter is phys- 
ically relevant only when |V$ N | < ao, where ao is the 
Milgrom acceleration p|, with ao ~ 10~ 10 ms _1 in order 
to fit galaxy observations without dark matter. The ex- 
tra field is governed by the non-linear Poisson equation: 



V • \pi (k \V(j>\ /a ) V0] = nGp 



(1) 
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where p is the matter density, k is a constant parameter 
and G is the underlying gravitational constant. On the 
left-hand-side, p(y) is an unknown function that must 
tend to 1 for j/ > 1 but behaves like y for y <C 1. We 
stress that this function is not the ratio of the Newto- 
nian and actual accelerations, which we will denote as 
p(x) with x = a/ao, a function more often employed 
by the astronomy community. Although for spherically 
symmetric cases the conversion between p.(x) and p(y) is 
straightforward, in more general situations complications 
may arise (see Pil). We will choose a particular form for 
fi(y) following Ref. [H (see Eq. M in Section |m] of the 
present paper), and likewise take k = 0.03 throughout 
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this article. We'll comment on other choices of // and 
their implications for our results, in the concluding Sec- 
tion of this paper. 

In regions where |V$ N | ^> ao, we have /i f=s 1 and 
therefore <j> yields accelerations that are k/Att times the 
Newtonian contribution. Hence by measuring the to- 
tal gravitational force we would measure the Newtonian 
gravitational constant to be (1 + k/Att)G and not G. 
However, given that here we take the value k = 0.03, 
we will ignore this negligible rescaling and assume G is 
the standard value measured by experiment. 

Analytical results exist for the above equation in the 
case of two point masses [l8| , But these are valid only in 
certain regimes. The "deep-MOND" solution of Ref. [H 
is valid close to saddle, while the "quasi- Newtonian" solu- 
tion is valid for large distances. By tackling the problem 
using numerical techniques we can present results accu- 
rate even in the gap between these two regimes. This 
is particularly important since this gap is in fact the re- 
gion of most practical importance for probes like LISA 
Pathfinder. 

Furthermore, since the Earth-Sun saddle point lies 
within the orbit of the Moon, we must concern ourselves 
with the effect of the Moon on the LISA Pathfinder mea- 
surements, and hence consider the more complex three- 
body problem. Our numerical approach enables us to do 
just that, and also to explore the possible advantages of 
probing the hitherto unexplored lunar saddle point. 

This article is laid out as follows. In Section [XT] we 
present the algorithm employed in our numerical code. 
Then in Section ITTll we test the code by comparing its re- 
sults with the analytical solutions previously developed 
for the two-body problem in [18j. Aware of the virtues 
and limitations of the code in Section IIVI we explore the 
solution for the full three-body problem, in a number 
of configurations, and considering both the Earth-Sun 
and Moon-Sun saddles. Specific recommendations for 
the LISA Pathfinder mission are made in Section [V] and 
more general considerations are included in the last Sec- 
tion. In Appendix [A] we present details of our numerical 
algorithm. 



II. NUMERICAL METHOD FOR SOLVING THE 
NON-RELATIVISTIC EQUATIONS 

In this section we present an overview of our numeri- 
cal approach for the solution of Eq. [T] (see Appendix [A] 
for full details). Most importantly this involves solving 
numerically only for the region immediately surrounding 
the saddle point, as shown in Fig. [T] This region contains 
no gravitational sources, so that Eq. ([1} becomes homo- 
geneous 1 . The gravitating sources make their presence 




FIG. 1: An illustration of the lattice location when calculating 
for the Earth-Sun saddle region. The calculation is performed 
for the a cuboid surrounding the saddle region, but not en- 
closing any of the gravitating bodies. Also shown is the type 
of non-uniform coordinates that we consider in order to boost 
the resolution in the central regions of the lattice, i.e. near 
the saddle point. For example, the y-resolution is a function 
of the y coordinate only. Not to scale. 

felt by the boundary conditions on our gridded volume. 
So long as the boundary is sufficiently far from the saddle 
point, so that 1 — fi <C 1, then to a very good approx- 
imation, (f> on the boundary is merely <f> — -p-<& , i.e. 
a rescaled version of the Newtonian potential. We also 
use the rescaled Newtonian solution inside the box as the 
initial configuration supplied to our numerical relaxation 
algorithm. 

It turns out that solving Eq. ([1} numerically in terms 
of <p is sub-optimal since it is only that appears in 
this equation. However, solving in terms of cither g = 
— V0 or u = — (Airfi/ K)Vcf) results in there being a further 
equation added to the system in order to specify the curl- 
free nature of g or u//i. Here we will use g, and so obtain: 

v • m = o (2) 

V x g = 0, (3) 

(written in terms of u, the y, appears instead in the sec- 
ond equation) . Fortunately we can bypass the extra work 
involved in dealing with the extra equation by choosing a 
relaxation algorithm for Eq. ^ that preserves the curl of 
g as it proceeds. Since we use a Newtonian initial config- 
uration, with zero curl, we obtain a solution to Eq. ([3]) at 
no additional computational cost. We represent the field 
g on a non-uniform lattice of discrete points, as shown 
in Fig. [T] Therefore, as we explain in Appendix [AJ in 
reality we obtain a solution to a finite difference version 
of Eq. @ . The algorithm preserves the finite difference 
version of V x g exactly (i.e. to floating-point accuracy 
in real life). 



This of course neglects the very small density present in inter- 
planetary space and any effects of the LISA Pathfinder measure- ment probe. 



3 



Interestin gly, this approach is the reverse of that used 



by Milgrom [21 



who solved for u using the Newtonian 
solution to obtain V ■ u = initially. The field u was 
then relaxed to a solution in which u/// has zero curl 
whilst not disturbing the divergence. However Milgrom 
considered only 2D problems, for which there is just one 
curl equation. In 3D there are three curl equations and 
it is beneficial to instead obtain the solution to the curl 
equations for free; hence our chosen algorithm. 



III. THE TWO-BODY PROBLEM 

We first consider a two-body setup, such as the Earth- 
Sun system in the absence of the Moon or any other per- 
turbations. As we shall show later, the Moon only per- 
turbs slightly the results for the Earth-Sun saddle point 
and hence this simpler problem provides great insight 
for the three-body situation to be explored in the next 
section. Furthermore there are analytical solutions for 
the two-body problem, so it presents a natural arena for 
checking our numerical algorithm. Conversely our nu- 
merical results will fix a poorly constrained parameter 
which appears in one of the analytical solutions. 

In Ref. [HI solutions in two different regimes were 
found for the MONDian two-body problem. The first 
is valid for regions far from the saddle point and was 
labelled the "quasi-Newtonian" solution; while the sec- 
ond is valid closer to the saddle and was dubbed the 
"deep-MOND" regime. Unfortunately from a practical 
perspective, it would be difficult to fly a spacecraft very 
close to the saddle point, where the deep-MOND solu- 
tion is accurate, due to navigational difficulties but also 
because, for example, the Earth-Sun saddle point moves 
as the Moon orbits (see next section) and is perturbed 
to a lesser extent by other solar system objects and by 
the Galaxy. In contrast, for large distances from the sad- 
dle, where the quasi-Newtonian solution is accurate, the 
MONDian effects are simply too small. Hence it is in fact 
intermediate distances, between the two regimes, which 
are of most practical interest and for which our numerical 
results have most value. 

In this section, unless otherwise specified, we will fol- 
low Ref. [l8[ and approximate the Newtonian accelera- 
tion (providing the MOND boundary conditions) by: 



(4) 



where A is a constant (related to masses of the Sun and 
the Earth and their separation; see [18|]), the problem is 
symmetric about the x-axis (connecting Sun and Earth) 
and the saddle is at (0,0,0). This is a valid approximation 
in the region near the saddle point and is used in deriving 
the analytical solutions. We shall use it for our numerical 
work in this Section (but not in the one that follows) in 
order to perform a fair comparison with the analytical 
results. 



We take the function /i(/c|V0|/ao) to have the form 
implied by: 



(•5) 



This is merely for illustration purposes: it was the form 
chosen in Ref. [l8[ since it eases analytical progress (with 
our choice of K = 0.03 also following that article). For 
the idealized case of spherically symmetric matter distri- 
butions this /i can be converted into the ratio of New- 
tonian to (total) MOND accelerations fi, which is often 
discussed in galactic astronomy (see eg. Ref. [1, [Io|, [H- 
[24[). However, galactic observations are of limited use in 
guiding our choice for \i since large uncertainties exist in 
the matter distribution within galaxies, which in practice 
will not be spherical, and a large extrapolation in accel- 
eration scale is required to obtain predictions for jx in the 
regime of interest here. The consideration of alternative 
fi functions and the conversion of a null result from a 
pass of LISA Pathfinder through a saddle region into a 
constraint on \x is left for future work [25j . 



A. Analytical solutions 

The analytical solutions are most easily expressed in 
terms of a dimensionless version of u = — (Air/j,/ K)V<f>, 
given by: 



U 



-/zV0. 



16n 2 aQ 47rao 
This may be also related to the acceleration via: 



(6) 



g = -W=^(l + ^) 1/4 ^, (7) 

for the chosen \x form. The field U must be divergence 
free in vacuum in order to satisfy Eq. (fT|) , while in order 
to keep 4> curl-free for this /i function, U must satisfy: 

A{l + U 2 )U 2 V x U + U x VC/ 2 =0. (8) 

The deep-MOND regime is characterized by U -C 1, 
when the (1 + U 2 ) term can be set to 1, whereas the quasi- 
Newtonian limit is characterized by U ^> 1. A transition 
region separates the two, located near the ellipsoid: 



2 , y 2 + z ' 2 



(9) 



where ro is the key length-scale of problem, given by 

167r 2 ao 



ro 



K 2 A 



(10) 



We have that ro ~ 381 km for the Earth-Sun saddle point 
(ignoring all other gravitating bodies). 
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1. Quasi- Newtonian solution 

In the quasi-Newtonian regime, well outside ellipsoid 
©, the solution can be conveniently decomposed as: 



U= ^N(^) + — B(V>) 
r r 



(11) 



where (r,ijj,9) is a spherical coordinate system centred on 
the saddle point with the two gravitating bodies located 
at ip = and tp = w. The first term is just the Newto- 
nian acceleration multiplied by K 2 /167r 2 ao, and the sec- 
ond "magnetic" term has finite curl and is sub-dominant. 
The angular profiles are given by: 



N(<0) = N r e r + N^ 
N r = -[l + 3cos(2^)] 

4 

= --sin(2-0). 



for the Newtonian component, and by 

B(V0 
B r 



B r (ip)e r + 5</>(V0 e V> 

2 7T 



B^p — 



5 + 3cos2?/> 3^3' 
tan" 1 ( V3 - 2 tan | ) + tan" 1 ( + 2 tan • 



(12) 
(13) 

(14) 

(15) 
(16) 



V^sm?/; 



7T cos -0 + 1 
3^3 sin-i/' 



(17) 



for the mag netic component (these formulae are all de- 
rived in |18j]). Note that it is only in this regime that the 
Newtonian-like and curl components can be split in this 
manner. 



2. Deep-MOND solution 

Inside the ellipsoid © the curl contribution becomes 
essential, and cannot be disentangled from the full field. 
A semi-analytical solution was found in Ref. [18J to be: 

U = c(^-)" _2 (F(^)e r + G(0)e v ,) (18) 

with a w 3.528, and angular profile: 

F(ip) = 0.2442 + 0.7246 cos(2t/>) + 0.0472 cos(4^) + . . . , 
G{ip) = -0.8334 sin(2V0- 0.0368 sin(4V>) + . . . , (19) 

which is almost identical to the Newtonian profile (F ss 
N r and G ~ N$). These formulae were derived in 

The normalization C is set by the boundary conditions 
and matching the two solutions suggests C ~ 1. Here 
we will determine C by comparison with our numerical 
results. 



B. Numerical results and comparison to analytical 
solutions 



We have first applied our numerical algorithm to the 
case where the Newtonian acceleration on the lattice 
boundary obeys the linear approximation to the two- 
body problem given by Eq. exactly as in the above 
analytical results. This will allow a fair comparison be- 
tween the two approaches. Considering the Earth-Sun 
saddle (tq ~ 381 km) and using a 257 3 lattice of phys- 
ical extent 10 000 km and central resolution ss 2.6 km, 
we find results for g as illustrated in Fig. [21 These can 
be seen to yield a good match to the two analytical so- 
lutions within their respective domains and provides the 
appropriate interpolation between them. 

Due to discretization asymmetries, the position of g = 
is not precisely at the central site of the lattice. The 
lowest value of |g| is instead found on a nearby site, but 
this is just a few kilometres away and represents good 
accuracy considering the 10 000 km box size. Before 
performing the comparison against the analytical results 
we first translate our numerical solution to take this offset 
into account, which can be important at very small r 
values. Without this, the gap between the thick solid and 
dashed lines in Fig. [5] would have been more noticeable. 

The numerical results enable us to determine the pre- 
viously poorly constrained C parameter appearing in the 
deep-MOND analytical solution. To determine this we 
first converted g into U via Eq. (JT]) and then deter- 
mined the ratio of the numerical results to the C = 1 
analytical values for all lattice sites within bounds of 
r/ro = 0.05 — > 0.5, chosen to be comfortably in the deep- 
MOND regime. We find a ratio of C = 0.839 ± 0.016 at 
each site, the central value from which we have used for 
all comparisons against the deep-MOND solution. 

However, g is not the key measurable quantity. We are 
instead interested in the tidal stresses, to which LISA 
Pathfinder is sensitive. When calculating the observ- 
able stresses we must subtract from g the unobservable 
rescaled Newtonian contribution (see discussion on the 
rescaling of G in the Introduction). Hence we introduce 
the notation: 



Si 



d 2 



d 2 <S> 



dxidxj Att dxidxj 



(20) 



for the observable MOND stress tensor, where Xi — x, y, 
or z. Under the linear approximation to the Newtonian 
acceleration field (Eq. [4j, the Newtonian stress tensor 
is just = A diag(l, —1/2, —1/2) and hence we must 
simply subtract a constant tensor from the raw MOND 
results. 

We illustrate the nature of this subtraction in Fig. [3] 
and Fig. 21 where it can be seen that the constant un- 
observable contribution dominates the full stress in the 
quasi-Newtonian regime (for diagonal elements of the 
tensor). As a result of this any inaccuracies in our nu- 
merical results become more important relative to the 
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FIG. 2: A comparison between the numerical and analytical 
results for components of g = — V</> in the two-body Earth- 
Sun case. Results are plotted as function of r for three pairs 
of ij) values: and 7r (top), 7r/4 and 37r/4 (middle) and ±7r 
(bottom), with symmetry relating the g-component values 
within in each pair up to the sign in the analytical case. In 
the numerical case this is also approximately true, except for 
very low r where the discretization asymmetry prevents it 
(this serves as an estimate of the discretization errors). Note 
that is zero analytically for ifi = and it/2 at all r, and 
that g r tends to zero for if) = 7r/4 at in the large r limit. We 
have used C = 0.839 for the deep-MOND analytical results 
in this figure. 
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FIG. 3: The full contribution to the two-body tidal stresses 
from the MOND field and the re-scaled Newtonian contribu- 
tion that must be subtracted from them in order to yield the 
observable component. The positive results are for dg x /dx 
while the negative ones are for dg y /dy. Results are shown 
for the line y — z — 0, for which the analytical results shown 
are symmetric to the interchange x — > —x. This symmetry 
is not quite realized in numerical case, as can be seen by the 
small differences between the dashed (x > 0) and solid lines 
(x < 0) for the numerical solution. The numerical stresses 
become unreliable for \x\ < 20 km. 



size of the signal, once the subtraction occurs. Further- 
more, when calculating the stress there are finite differ- 
encing errors, which become increasingly important at 
small |a;|. This is evidenced by the differences between 
the dashed (x < 0) and solid (x > 0) lines for the nu- 
merical results, and their difference relative to the deep- 
MOND analytical result. Hence neither regime is trivial 
computationally. However, as already noted, for LISA 
Pathfinder we are fortunately interested in the interme- 
diate regime where the reliability of results is greatest. 

When a real spacecraft performs a saddle fly-by, it will 
not pass precisely through the saddle point. We dis- 
cuss the likely impact parameter and trajectory for LISA 
Pathfinder in Sec. |Vj but for illustrative purposes we con- 
sider here the form of S yy and the numerical uncertainties 
in it for a trajectory along the line y — 100 km. While the 
normal Newtonian signal is much larger than that from 
</>, the former is simply constant (though it would slowly 
drift without the linear approximation used here). The 
MOND signal, on the other hand, provides a distinctive 
variation as the probe passes by the saddle. The form is 
stable against numerical errors, although minor asymme- 
tries can be seen near the point of closest passage. These 
are slightly more significant when the lattice spacing is 
increased by enlarging the box side to 20 000 km (from 
the 10 000 km used thus far). However, the larger box 
provides greater accuracy in the less important tails of 
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FIG. 4: The observable contribution to the dg y /dy tidal stress 
from the MOND field for line y = z — 0. It can be seen that 
the the numerical results perform well even after the subtrac- 
tion of the unobservable re-scaled Newtonian contribution, 
which heavily dominates the stresses in the quasi-Newtonian 
regime. However, near the edge of the simulation volume the 
stress is attenuated due to the smooth transition of 4> to meet 
the rescaled-Newtonian boundary conditions. 



1 





2 


x10" 14 


















1 


CVJ 
CO 








10" 


2 










CM 




numerical (10 000 km) 










3 


numerical (20 000 km) 










analytical -500 


500 






Standard Newtonian / 


i x / km 






1 0" 


4 










10" 


5 










10" 


6 


^^^^ / 
~*-r-T, ^T. . ... i ........... , 









-8000 -6000 -4000 -2000 2000 4000 6000 8000 
x / km 



FIG. 5: The two-body MOND stress signal S yy along the line 
y = 100 km for two simulation box sizes (10 000 km and 
20 000 km) with fixed lattice size 257 3 . The zoomed region 
with linear y-axis shows the minor differences between these 
two sets of results near the point on this trajectory that is 
closest to the saddle, in addition to minor asymmetries which 
are due to discretization errors and more significant for the 
larger box. The log-scale graph shows the tails of the signal, 
which are improved in the larger box. 



the signal (see Fig. [5]). 

To close this Section we present results using the full 
two-body Newtonian acceleration field for the Earth-Sun 
system, rather than the linear approximation (Eq. (Q}). 
As can be seen in Fig. [51 the effect of the approximation 
is insignificant. Fig. [6] also shows the change in the form 
of the signal for three different impact parameters: y — 
25, 100 and 400 km. We see that the signal broadens 
with increasing y, but more importantly, its amplitude 
decreases. 



Sun-Earth- Moon system presents only two saddle points: 
the Sun-Earth saddle, which is more or less as the two- 
body analysis suggests, and a single saddle point close to 
the Moon, which is intrinsically a three-body problem. 
But even for the latter, we shall empirically find a recipe 
for adapting the two-body results to the real situation. 



A. The Sun-Earth saddle 



IV. THE THREE-BODY PROBLEM 

In this section we add a third body to the problem and 
consider the perturbing effect of the Moon on the Earth- 
Sun saddle. We also present results for the lunar saddle 
point (which is intrinsically a three-body problem, as we 
shall see). 

We must first comment that the "realistic Solar sys- 
tem" , as described in Ref . [l8[ , is in fact not realistic at 
all, and numerical work is needed to provide even quali- 
tatively correct directions for a space mission. The argu- 
ment in Ref. 18[ regarding the Earth-Moon saddle, for 
example, is entirely flawed because of its use of the lin- 
ear approximation for the Newtonian field for the Earth- 
Moon system while ignoring the gravitational effects of 
Sun (cf. Eq. (72) in that paper). As it happens, the 
perturbation induced by the Sun is too large for this ap- 
proximation to be valid and in fact there is no true Earth- 
Moon saddle point. Instead as explained in Fig. the 



For the Sun-Earth saddle, the effect of the third body 
is not drastic, because the lunar mass is approximately 
l/81th that of Earth and even at new Moon, the Earth 
is only about twice as far away from the saddle point. 

We will firstly consider the location of the saddle, 
which is of great practical importance when planning a 
spacecraft fly-by. It must be considered that the length 
scale on which the MOND signal is significant, ro, is ap- 
proximately 381 km (see Eq. [9]), so the saddle has to be 
located to a precision better than this. The saddle point 
is approximately 258 800 km from Earth (around two- 
thirds to the lunar orbital radius), and therefore even a 
small perturbation is potentially important. As shown 
in Fig. [51 at full Moon the saddle shifts about 250 km 
towards the Sun. The magnitude of the shift remains 
at approximately this level until the crescent phase ap- 
proaches, when the Moon begins to approach the saddle. 
The saddle then starts to experience a large perturba- 
tion and, in the 3 day period surrounding new Moon, 
it quickly moves over the left half of the quasi-ellipse 
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FIG. 6: The MOND stress signal S yy along the lines y — 
25, 100 and 400 km (top to bottom) for the full two-body 
Earth-Sun case (thick black dashed) and when using the linear 
approximation to the Newtonian acceleration (black solid). In 
the upper panel the rescaled Newtonian stress is shown (grey) 
for the y — 25 km case. 
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indicated in the figure. Note that the effect near the 
new Moon is significantly enhanced if the Moon is at 
its perigee during this period. While the perturbation is 
always less than 10 000 km, and is less than 1 000 km 
for most of the lunar cycle, this cannot be ignored if a 
spacecraft is to be navigated through the saddle region. 

While the perturbation caused by the Moon generally 
breaks the axial symmetry present in the two-body dis- 
cussion of the previous chapter, during a New or Full 
Moon the symmetry remains unbroken. Hence we may 
consider these two extreme cases via the same linear ap- 
proximation to the Newtonian acceleration, but with dif- 
fering A values (cf. Eq. (Q}). We can then make the 
following scaling argument. The conditions on the field 
U are that it must be divergence-free, satisfy Eq. (JSJ and 
match the boundary conditions. The first two of these 
are unaffected by a rescaling of the spatial coordinates 
r — > ar while the boundary conditions are such that: 



FIG. 7: Maps of the Newtonian |V$ N | for different con- 
figurations of 3 point masses, showing differing numbers of 
zero solutions. Black regions denote low |V$ N | values while 
white denotes high values. The upper-left map is for the case 
of equal masses arranged in an equilateral formation, which 
yields 3 saddle points plus a central maximum in $ N . The 
second frame shows the same number of solutions, but is for 
a configuration much closer to the Sun-Earth-Moon system, 
albeit that the Earth's mass is distributed evenly between 
the lighter two bodies and the heavier body is slightly lighter 
than the Sun. Mass is added to the heavier body in the third 
frame, causing the lower saddle and central maximum to ap- 
proach each other. Adding yet further mass causes them to 
meet and disappear, leaving only two saddle points by the 
time the solar mass is reached, as in the fourth frame. Giv- 
ing the lightest two bodies masses corresponding to the Earth 
and the Moon, then maintains this number of saddle points 
but, as can be seen in the final frame, the lunar saddle point 
is surrounded by only a small region of low acceleration and 
it is also significantly tilted towards Earth. 



g 



A-kclp, 



x . 1 y . 1 z , 
— l j k 

r 2 r 2 r 



(21) 



in the linear approximation of Eq. (j4]). Hence in scaled 
coordinates r/ro there is a single solution for U valid for 



327 km 2 . Hence the Newtonian stresses will be slightly 
higher at new Moon than in the two-body case, and very 
slightly lower in the Full Moon case. Considering these 
changes at a fixed r/ro value, the MOND g will be un- 



all A. Therefore while the form of S« stresses is unaf- affected ' therefore 5 u at fixed r/r increases in inverse 

fected by the change in r„, the spatial extent of the signal Proportion to r„. However, note that at new Moon, the 
is proportional to r$ and the stress magnitude scales in- 
versely with ro. (This is a general argument to be used 
later in the Moon saddle analysis). 

i , r, „ , , , . . . . . These figures assume the Moon-Earth separation is equal to the 

At .bull Moon the change m rn is minor, increasing . . . c , ... r , . . 

° u ' ° semi-major axis tor the lunar orbit; tor the lunar perigee this 

by about 1 km from 381 km. However at new Moon, value should be decreased by 22 km and for the apogee it should 

the change is more significant, with ro decreasing to be increased by 15 km. 
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Ax / km 

FIG. 8: The perturbation in the position of the Earth-Sun 
saddle due to the presence of the Moon when the lunar orbital 
radius is fixed at its semi-major axis (thick), with addition 
plots for apogee (outer) and perigee (inner) to highlight the 
effect of orbital ellipticity. Negative Ax values denote that the 
saddle moves closer to Earth, which would lie to the left in this 
figure, while the Sun would lie to the right. The saddle is only 
slightly perturbed at Full Moon, moving slightly towards the 
Sun, but as new Moon approaches the perturbation becomes 
large and saddle moves quickly, covering the left half of the 
quasi-ellipse in approximately 3 days, with new Moon being 
the mid-point of this period. The perturbation at new Moon 
is enhanced if new Moon coincides with the lunar perigee but 
attenuated if at apogee. 



linear approximation is not so robust because the gravity 
field of the Moon varies on a smaller length scale than 
that of the Earth. 

To consider an arbitrary Moon position, and without 
the approximation to the Newtonian acceleration field, 
we require numerical methods. The results in Fig. [51 
which shows three different lunar phases, indicate that 
the effect of the Moon on the MOND signal near the 
Earth-Sun saddle point is minor. At Full Moon the signal 
is essentially identical to the two-body case, since the 
Moon is more than twice as far from the saddle point as 
Earth. The saddle is displaced furthest away from the 
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FIG. 9: The MOND stress signal S vv along the lines y = 25, 
fOO and 400 km (top to bottom) for the three-body Earth-Sun 
case for different lunar phases: new Moon (thick, black, solid), 
Full Moon (thick, black, dashed) and when the Moon appears 
18° away from the Sun towards positive y (thin, black, solid). 
Also shown in the y = 25 km case are the rescaled Newtonian 
stresses (grey). 

Earth-Sun line when the angle between the Moon and 
the Sun (measured from Earth) is approximately 18°, 
and even then the line in the figure that corresponds to 
this asymmetric case still has the same basic form as the 
two-body signal. Some x —> —x asymmetry is instilled, 
but the overall conclusion is that the MOND signal for 
the Earth-Sun saddle is largely unaffected by the phase 
of the Moon. 



B. The Lunar saddle 

The lunar saddle can potentially provide significantly 
larger stress signals than the Earth-Sun saddle and is 
therefore of potential interest experimentally. However, 
the larger stress signal implies that the MOND region is 
smaller and so the challenge of navigating a space probe 
through the saddle region becomes more difficult. 

1 . Newtonian gravity field 

We must first consider the lunar saddle from a purely 
Newtonian perspective, since the Newtonian physics sets 
the saddle location and the boundary conditions for our 
MOND calculation. It is useful to think of the lunar 
saddle as a heavily perturbed Moon-Sun saddle, since 
the Earth is the least important of the three bodies in 
the saddle vicinity. 

Although the gravitational pull of the Sun is effectively 
constant on the scale of the lunar orbit, the Moon is ap- 
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proximately l/81th of the mass of Earth, and therefore 
the distance from the unperturbed Moon-Sun saddle to 
the Moon should be around l/9th of the distance from 
the Earth-Sun saddle to the Earth. Since tidal stresses 
follow an inverse-cube law then the Sun makes only a mi- 
nor contribution to the tidal stresses at either saddle, but 
the Moon causes stresses at the unperturbed Moon-Sun 
saddle around 9 times stronger than the tidal stresses 
caused by the Earth on the Earth-Sun saddle. This im- 
plies (cf. Eq. ITOl) that the value of r$ for the lunar sad- 
dle should be around 9 times smaller than that for the 
Moon-Sun saddle, before third-body perturbations are 
considered. 

However the perturbing effect of the Earth on the lu- 
nar saddle is very large, as can be seen in Figures [TOl and 
fTTI At new Moon the Earth helps the Moon overcome 
the pull from the Sun and hence the saddle is further 
from the Moon than on average and the lunar stress con- 
tribution is smaller, yielding a larger saddle region. At 
Full Moon, the Earth acts with the Sun and therefore 
the saddle is closer to the Moon, the stresses are greater 
and the saddle region becomes very small. Additionally, 
the Earth shifts the lunar saddle significantly away from 
the Moon-Sun line, by up to ps 30° (measured at the 
Moon), although of course it lies on this line at full and 
new Moon (ignoring the possible small displacement of 
the Moon off the Ecliptic) . The lunar saddle is therefore 
intrinsically a three-body problem. 

Importantly, since the Moon dominates the tidal 
stresses near the saddle, the form of the Newtonian accel- 
eration in the saddle region is still approximately that of 
the linear approximation form, just as in the Earth-Sun 
case where the Earth dominates the stress. However, now 
the effective "axis of symmetry" points from the Moon 
to the saddle location. At once this suggests a recipe for 
adapting the two-body MONDian solution to the lunar 
saddle. 



2. MOND stress signal 

Given the observation that the Newtonian acceleration 
field near the lunar saddle is roughly linear (as in Eq. [4j , 
it might be expected that the MOND solution would be 
very similar to the two-body case, except for re-scalings 
with tq. This is indeed seen in the numerical results 
displayed in Fig. [121 where the new Moon case is indis- 
tinguishable from the previous Earth-Sun results and the 
even the non-axial quarter Moon case is only slightly dif- 
ferent, once appropriate scalings with ro are introduced. 
For the purposes of this figure we have introduced a new 
coordinate system (x',y',z), which is a rotated version 
of (x, y, z) such that the x'-axis joins the Moon and the 
saddle point, which still lies at (0,0,0). We have then 
chosen y' ss 0.26ro and 1.05ro lines for which to plot the 
MOND stress signals, where ro ~ 81 km at new Moon 
and ro ~ 38 km at quarter Moon. The two panels in this 
figure are hence directly comparable with the lower two 




FIG. 10: Upper panel: The distance of the lunar saddle from 
the Moon as a function of the lunar phase, where 9 m00 n = 
denotes new Moon and # moon = 180° denotes Full Moon. 
Middle panel: The effective ro value for the lunar saddle as a 
function of lunar phase. Lower panel: The deflection angle of 
the saddle away from the Moon-Sun line (towards Earth). 
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FIG. 11: Four maps of the Newtonian acceleration strength 
|g N | (as before dark means small |g N |, and so a saddle, in- 
dicated by a X wherever it makes it clearer) in the vicin- 
ity of Earth for Sun-Earth-Moon system (in plane 2 = 0). 
Clockwise from top left the angle between the Moon and the 
Sun, measured at Earth are 0° (new Moon), 18°, 90° (First 
Quarter) and 180° (Full). The Earth (0,0) and the Moon are 
shown by crosses. White corresponds to values greater than 
80 x 10 6 ao, black indicates values less than 25 x 10 6 ao, with 
shades of gray indicating intermediate values. Note that the 
Earth perturbs the Moon-Sun saddle region heavily, both in 
orientation and size. 
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FIG. 12: Stress signals for the lunar saddle, compared to the 
results from the Earth-Sun saddle (in the linear approxima- 
tion). Lunar results are expressed in a new coordinate system 
(x',y',z), a rotated version of (x,y,z) such that the x'-axis 
is the line joining the Moon and the saddle, with y ~ 0.26ro 
(top) and 1.05ro (bottom). In the Earth-Sun case, the re- 
sults are for y — 100 and y — 400 km, which correspond to 
the 0.26ro and 1.05ro, making the figure directly comparable 
with the lower two panels of Fig. [6] At New moon ro ~ 81 km 
while at Quarter Moon ro ~ 38 km. 



35 cm apart. These test masses will follow geodetic mo- 
tion to exceptional precision due to protection from solar 
radiation pressure and other unwanted non-gravitational 
influences, while the spacecraft itself is designed to min- 
imize its own gravitational impact. An interferometer 
will study the relative motion of the test masses such 
that LPF will be an excellent instrument for measuring 
tidal stresses, such as those predicted in the preceding 
two sections. The mission is presently in the "implemen- 
tation" phase with its launch expected in 2012. 

The nominal plan for the LPF mission is for it to en- 
ter a large-amplitude Lissajous orbit (of near halo-orbit 
dimensions) around the LI Lagrange point of the Earth- 
Sun system, at a distance of ~ 1.5 x 10 6 km from Earth. 
Once the planned mission at LI is complete, it can in 
principle be extended, by breaking from the Lissajous 
orbit and passing close to the Earth-Sun saddle point, 
or alternatively near the lunar saddle. In this section we 
compare the sensitivity of LPF to the possible MOND 
signal, given its navigation and measurement capabili- 
ties. 



panels in Fig. [51 However, we have not included results 
equivalent to the upper panel there (y' « 0.066r ) since 
at Quarter Moon this amounts to y ~ 2.5 km which we 
believe is to too small an impact parameter to consider 
for LISA Pathfinder (see next section). 

Hence the form of the lunar saddle MONDian signal 
is essentially the same as for the Earth-Sun saddle; yet 
it offers larger MOND and Newtonian stresses, but over 
a smaller region. In practical terms, however, not much 
benefit might arise if the navigational constraints dictate 
a certain minimum likely impact parameter. The MOND 
signal varies roughly as 1/r, so for this impact parameter 
the likely stress level is largely independent of r . 

V. TESTING MOND WITH LISA PATHFINDER 

We now explain more specifically how the results found 
in this paper translate into a concrete signal, given the 
navigational constraints and instrument noise properties 
of LISA Pathfinder. We consider one example here and 
in a companion paper (26| perform a more comprehensive 
study of the issues highlighted. 

A. Overview of LISA Pathfinder Project 

The LISA Pathfinder (LPF) project [H is an Euro- 
pean Space Agency (ESA) mission designed to test the 
technology intended to be employed in LISA (Laser In- 
terferometer Space Antenna), a proposed gravity wave 
observatory. While LISA itself is planned to consist of 
three spacecraft in a triangular configuration of side- 
length 5 x 10 6 km, the Pathfinder mission will consist 
of two test masses within a single spacecraft and just 



B. Spacecraft navigation 

The objective of the LPF mission does not require the 
spacecraft to have the large thrusters of an interplane- 
tary probe, since the journey out to LI is made possible 
by a propulsion module from which LPF then separates. 
Instead LPF will carry only low-thrust micro-propulsion 
systems and the ability to navigate the probe near a sad- 
dle point is not something that can be taken for granted. 

Fortunately, the Lissajous orbit is unstable and hence 
LPF can be ejected from the LI region using a fairly 
small change in its velocity. For example, a 30-day burn 
can yield a dv of ~ 1 ms" 1 and put LPF on a trajec- 
tory that brings it close to Earth in a reasonable time- 
frame, of order one year. The precise orbit depends on 
the timing of the burn relative to the phase of the Lis- 
sajous orbit, the burn magnitude, subsequent correction 
manoeuvres and any close lunar fly-bys. Without the 
latter two, approaches near the Earth-Sun saddle in ap- 
proximately 1 yr tend to have quite large impact param- 
eters ~ 10 000 km but closer approaches < 1 000 km 
can be achieved using a lunar fly-by. Further correction 
manoeuvres can aid targeting the saddle region and the 
likely impact parameter is ultimately determined by the 
ability to track the spacecraft and apply appropriate dv 
manoeuvres. Due to Earth's gravity the spacecraft will 
be cruising through the saddles with a speed of the or- 
der of ~ 1 kms . In what follows we will take it that 
LPF can be directed through a saddle region with veloc- 
ity 1 kms -1 and an impact parameter ~ 50 km. A more 
complete study of these orbits is presented in [26| . 
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C. Sensitivity to MOND stress signal 

Firstly, it should be noted that the design of LPF lim- 
its it to being sensitive only to one of the diagonal el- 
ements of the stress tensor, eg. S yy in our notation. 
Additionally, in order to yield a stable radiation pres- 
sure it is desirable to keep fixed the orientation of the 
(roughly cuboid) probe relative to the Sun, with the side 
on which the solar panel is mounted facing directly at the 
Sun. Given the vast distance of the spacecraft from the 
Sun, the solar panel is therefore effectively aligned nor- 
mal to the x-axis and the spacecraft design then implies 
that S xx cannot be measured. This leaves the measure- 
able stress as cos 2 (a)S yy + sm 2 (a)S zz , where a specifies 
the orientation of the test masses relative to the y-axis. 

Estimates of the sensitivity of LPF as a gradiometer 
are normally expressed via the power spectral density of 
the stress signal. This is defined as: 



P(f) 



2 
T 



-T/2 



-T/2 



(22) 



where / is the frequency, t is the time and T is the in- 
tegration period. When integrated over / from to in- 
finity this yields the mean "power" in the signal, i.e. the 
time-average of Sfj , with the factor of 2 present because 
negative frequencies are folded in with the positive ones. 
For the present article we assume that LPF will attain a 
sensitivity in the square-root of this quantity, the ampli- 



tude spectral density, of 1.5 x 10~ 14 s~ 2 /vHz between 1 
and 10 mHz [20(. Beyond this range we will take the sen- 
sitivity to degrade as 1// at lower frequencies and as / 2 
at higher frequencies. This yields a good approximation 
to published expectations [20| . 

We simulate gradiometer noise under the assumption 
of Gaussianity (independent Fourier mode phases) and 
the above spectrum. We add this to the theoretical stress 
signal, which is the combination of the standard Newto- 
nian and MOND components, and then apply a simple 
cleaning algorithm to extract a likely measured MOND 
signal. Specifically, we approximately remove the New- 
tonian signal by performing a quadratic fit to the data. 
In the linear approximation, the Newtonian stress sig- 
nal is just a constant while the gravitational acceleration 
near the saddle is so small that we may assume the probe 
moves with constant velocity and therefore this fit takes 
into account two further orders, beyond linear, in the 
spatial Taylor expansion of the Newtonian acceleration. 
We then apply a top-hat band-pass filter to remove all 
frequencies outside the range 0.1 to 10 mHz, i.e. keeping 
those of the approximate order as the MOND signal. 

The results in Fig. [T3] assume a spacecraft velocity of 
1 kms -1 through the Earth-Sun saddle region, passing 
on a trajectory with y = 50 km and z = 0, during a full 
Moon. While the velocity of the probe is unlikely to be 
aligned with the x-axis in this manner, the details of the 
trajectory orientation do not appear to be greatly impor- 
tant. A significant signal is clearly seen in both real-time 
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FIG. 13: Theoretical and recovered stress signals S yy and S zz 
for a y — 50 km, x = 1 kms -1 trajectory, including a real- 
ization of gradiometer noise and applying a simple cleaning 
algorithm (see text). These results are for the three-body 
Earth-Sun saddle with the lunar phase being full Moon. 



and Fourier domains: Figs. [T3] and Q3J respectively. If 
the actual trajectory saw LPF cross the saddle region 
more quickly, then the MOND signal would be shifted to 
higher frequencies where the signal to noise is even bet- 
ter. An increase in the impact parameter r m ; n reduces 
the signal roughly as l/r m ; n for r m ; n ~ tq, but as 1/r: 



2 

min 



for r min > r . 



VI. CONCLUSIONS 



We solved the non-relativistic limit of TeVeS for two 
and three-body saddle points using a particular [i func- 
tion used previously in the literature. We predicted an 
anomalous stress signal that is potentially detectable by 
future space-probes. We found that past analytical calcu- 
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FIG. 14: The amplitude spectral density (ASD) for the S zz 
MOND signal for a y = 50 km, x — 1 kms -1 trajectory, 
compared to the simple noise model employed here. Note 
that while the noise ASD is independent of the integration 
time, there being a single MOND event means that its ASD 
lessens with increasing integration time, where for the present 
plot we conservatively integrate over 2 x 10 4 s. 



lations for the two-body case are not only in agreement 
with our numerical results, at least in their regimes of 
applicability, but provide a good insight into three-body 
problems of the Sun-Earth-Moon type, where the New- 
tonian stresses at the saddle points are dominated by a 
single body. Our results indicate that the most significant 
effect of the Moon on the Earth-Sun saddle is to change 
its location, while the MOND anomalous stress signal is 
not greatly affected, and is stable to changes in the phase 
of Moon. Furthermore, we have demonstrated that the 
single lunar saddle point near the Moon yields a stress 
signal of similar form but differs in stress magnitude and 
spatial extent. However, for an impact parameter likely 
from a spacecraft flyby, we note that these two factors 
effectively cancel (at least in some regime). 

Finally, we have demonstrated that a possible exten- 
sion to the LISA Pathfinder mission, targeting the sad- 
dles, would be sensitive to the anomalous stress signal 
that we have calculated. This is in part due to the dis- 
tinct form of the signal close to a saddle point, which 
cannot be confused with the (much smoother) Newtonian 
signal (or with any other signal derived from the param- 
eterized post-Newtonian framework). It is also due to 
the likely speed of the probe through the saddle region 
yielding signal variations at frequencies that minimize 
the instrument noise. While this conclusion assumes the 
ability to pass within ~ 100 km of a saddle point, we 
believe that this is achievable even with the low thrust 
propulsion system incorporated in the LISA Pathfinder 
design. In future work we hope to implement our studies 
as part of a concerted effort to explore realistic trajecto- 
ries for LISA Pathfinder. 



We conclude with a general remark on MONDian the- 
ories which we hope will clarify a number of issues with 
the proposed test. For all MONDian theories with a rela- 
tivistic formulation (and not only TeVeS, as described in 
the introduction), MONDian non-relativistic effects are 
due to an extra field, 4>, as opposed to the gravitational 
field itself, i.e. the metric field, with 3"^ associated with 
goo. Indeed it's very hard to covariantly modify a gravity 
theory so that in the non-relativistic limit g 00 is ruled by 
a non-linear Poisson equation. In the Newtonian limit 
(when fj, w 1 in Eqn.([T])) the field <f> mimics the New- 
tonian field and so renormalizes the gravitational con- 
stant (G — > G(l + k/(47t)) as discussed in the literature 
(e.g. pl). To satisfy constraints, <j) must be sub-dominant 
in this regime, and this is enforced by the parameter k 
appearing in its Poisson equation ([1]). 

But this simple fact implies that one must trigger 
MONDian behavior in <j> at Newtonian accelerations a^v 
much higher than eto. Only thus can </>'s relative impor- 
tance start to increase with decreasing a at, so that by 
the time a at ~ »o the field <f> is not only MONDian but 
dominates (&n, as required by astronomical applications. 
It can be easily computed that if fj, turns from 1 to a 
single power-law, then MONDian behavior in <f> should 
be triggered when the Newtonian acceleration drops 
below (47r/K) 2 a ~ 1.75 x 10 5 a ~ 10~ 5 ms~ 2 . And in- 
deed this is the rough Newtonian acceleration at r$. It 
is therefore important to realize that for LPF realistic 
impact parameters we would probe the regime where <fi 
has gone fully MONDian, but hasn't yet dominated <J>at. 
The MOND signal can be detected above the Newtonian 
one because it has a distinctive spatial variation whereas 
$jv is just a DC component. 

Bearing this in mind a number of issues may be clari- 
fied. One concerns the self-gravity of LPF. This is only 
balanced at the level of a at ~ 10~ 9 ms~ 2 , so one might 
think that a test of MOND, in the regime a at ~ a>o, would 
run against the wall of self-gravity. In fact we are testing 
much higher Newtonian accelerations; for instance, for an 
impact parameter of 40 Km we have ~ 10~ 6 ms~ 2 . 
We'd need to approach the saddle much closer than about 
400 meters before self-gravity became an issue (and the 
spacecraft itself had to be included in the computation 
of the location of the saddle) . 

This should also clear the matter of the generality of 
the predictions made. There are several /i on offer for 
astrophysical purposes. They all have the same rough 
behavior in the regime under study, where <f> is fully 
MONDian but hasn't yet dominated $at. A negative 
result from LPF would rule out virtually all proposed ji. 
It would only fail to rule out very contrived /j,(y) (never 
suggested in the literature), for which two power laws are 
employed in /i(y): a very steep one from y = 1 (felt at r ) 
to the point where a at = ao; and the usual linear power- 
law for ajv < clq- extremely contrived. We are currently 
writing a follow up paper expanding on this matter, 

The high sensitivity of LISA Pathfinder enables it to 
potentially detect the anomalous MOND stress when the 
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deviations from Newtonian dynamics are tiny, in con- 
trast to the conditions present in outer reaches of galax- 
ies 0, HI]. On the other hand, constraints are yielded 
by planetary orbits [l5j only at much stronger accel- 
erations than LISA Pathfinder would probe. A saddle 
point fly-by offers the chance to study MOND where we 
have an exquisite knowledge of the mass distribution and 
the ability to make full calculations of the MOND accel- 
eration field, without assuming simplifying symmetries. 
Furthermore, saddle point exploration would enable an 
investigation of gravity at very low accelerations in the 
close vicinity of Earth, which can be reached in a short 
time-frame and without the increased mass distribution 
uncertainties affecting probes sent to very large distances 
from the Sun. A positive detection would, of course, vin- 
dicate the MOND paradigm, while a null result would 
provide clean and reliable constraints upon it. 

In closing we note that several refinements to our 
numerical calculations should be trivial. Including the 
quadrupole and higher multipole moments of the sources 
or the effect of Jupiter and other Solar system objects (or 
even the galactic field) is straightforward, as they only af- 
fect the boundary conditions for our problem. However 
these effects are expected to be small. The shift of the 
saddle location is of course non- negligible, but can be eas- 
ily computed with effects already studied in [lis}. Once 
the new location is taken into account the change in the 
detailed MOND tidal stresses can be computed just by 
changing the boundary conditions of our code. However, 
if the effect of the Moon on the Earth-Sun predictions 
for realistic impact parameters is already so small, these 
effects can be expected to be completely negligible. 
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Appendix A: Details of numerical method 

Here we briefly describe our numerical algorithm. This 
is based upon representing g on the sites of a non- uniform 
lattice of the form shown in Fig. [1] such that we obtain 
greater spatial resolution near the saddle point. A relax- 
ation algorithm then cycles over each lattice site x and 
changes the values of g at x and the neighboring sites 
so that the divergence equation Eq. © is solved locally, 
that is: 



A, 



E 



o. 



(Al) 



We then move to the next site and change the field val- 
ues so that the condition is valid there, using the newly 
updated values as we proceed. However, when enforcing 
the above condition at these subsequent sites, the value 
of Z? x at the first site will be slightly changed. We there- 
fore require many cycles over the whole lattice before the 
above condition is closely matched globally. 

We stress that we use a different discretization proce- 
dure to Ref. [2l|, in that we define all three components 
of g as well as /i at the same location, that is: 



/i x = /i(Kg x /a ). 



(A2) 



This is in contrast to the more elaborate scheme used in 
Ref. [2l[, in which u 3 is defined on the j-links between 
the lattice sites and then [i is defined at the centers of 
the grid squares in the 2D case considered. Under that 
scheme the calculation of u l j ' [i at any position in a 3D 
calculation would require knowledge of 33 values of u 3 , 
whereas to determine /u x g x here, we require merely the 
knowledge of the three components of g x . 

The presence of the non-linear function /i is a no- 
table complication. Furthermore its form is not uniquely 
known and therefore we desire to be able to solve for an 
arbitrary function. The algorithm of Ref. [2lJ involves 
solving the curl equation for u locally, while keeping the 
\i values fixed at their old values. Once the [i values 
are updated using the new u, the curl equation is no 
longer matched, and this slows the convergence of the 
algorithm. In contrast, here we proceed by solving the 
divergence equation to first order in <5g and <5/i, where S 
denotes the change from one step to the next. Then, as 
the system converges to the solution, the terms of order 
6 2 become negligible very rapidly. 

Crucial in the above is to ensure that in updating the 
field configuration the curl of g remains zero. We define 
the our preserved discrete curl as: 
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where r x is the position vector site x. In order to preserve 
this discrete curl at each step of the relaxation we change 
the fields according to: 



5gi = 



^x-j 



(A4) 
(A5) 



where the value of C x is chosen, as below, to yield 
Eq. (|Alj) to first order. While the first equation acts 
on all three components of g, the second acts on only 
one component per site. 

To determined C x , consider that to first order in 5g 
and 5[i: 



E 

j 



Aix<5.9 x + gi8^ - Hx-jSg 3 - - ,g x :<5^x-j 



.(A6) 



x-J 
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Writing 6/j, in terms of 5g 2 and substituting the above 
values for <5g then yields: 



6D» 



E 

3 



2E 



AV A J 1 d 5 2 



Mx-j 
{MY 

where we have used compact notation such that: 




A J = 



yj - 



x-J' 



(A7) 
(A8) 



(A9) 
(A10) 



We can then insert the fi derivatives for our chosen /i 
function. Finally, since we want D x to be zero after each 
change, we set: 



6D X = -A, 



(All) 



and hence know the C x required to obtain D x = 0, at 
least to first order. (If had not changed during this step, 
the above procedure would have set D x to be exactly 
zero, i.e. to all orders in Sg.) 

Note that in addition to preserving the curl, Eqs. (IA4|) 
and (|A5j) also exactly preserve the change in </> measured 



across the lattice as: 



A<t> = 



(A12) 



where x % and x k are fixed during the summation. 

In practice, cycling over the lattice and solving the 
discrete equation locally does not lead to rapid enough 
convergence to the solution. As D x is (approximately) 
zeroed at later sites, D x at earlier sites is moved slightly 
away from its desired value and a large number of it- 
erations of this procedure are required before each g x 
has converged to a good approximation. We there- 
fore preempt the changes to the field that will occur at 
other points in the cycle using the fact that these very 
changes are (largely) responsible for D x being non-zero. 
This is achieved by a method known as successive over- 
relaxation (SOR, e.g. [13]) in which Sgi — > XSgi , where A 
is the over-relaxation parameter and is larger than unity. 
We begin with A = 1 and increase it once the system has 
begun to settle down, since high values of A can initially 
result in the RMS value of D x increasing, contrary to our 
goal. 

We coded the algorithm outlined in this Appendix us- 
ing the LATfield library [HI]. 
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